CRB Variance

Disclaimer

This work is very preliminary as I get back into the coding swing of things. Data wrangling and figure generation will be done via R, but the rest of the project will be done using good ol’ microsoft products. This is just an entry point into data crunching and should by no means be considered a final product.

Steamboat Tower SNOTEL

Raw data download

Steamboat Tower (825) SNOTEL data

# The fucntion snotel_download evidently downloads everything in metric units. 

snotel_825_tower <- snotel_download(site_id = 825, path = tempdir('../data'), internal = TRUE)

## Error in wdman::phantomjs(verbose = FALSE) : 
##   PhantomJS signals port = 4567 is already in use.

write.csv(snotel_825_tower,"C:/Users/13074/Documents/ESS580/thesis_project/thesis_research/data_raw/snotel_825.csv", row.names = FALSE) #write in the raw data

Data Cleaning

head(snotel_825_tower) # check the date, usually a character.  

##   network state site_name               description      start        end
## 1    SNTL    CO    tower  Fish Creek (140500010407) 1978-10-01 2022-06-15
## 2    SNTL    CO    tower  Fish Creek (140500010407) 1978-10-01 2022-06-15
## 3    SNTL    CO    tower  Fish Creek (140500010407) 1978-10-01 2022-06-15
## 4    SNTL    CO    tower  Fish Creek (140500010407) 1978-10-01 2022-06-15
## 5    SNTL    CO    tower  Fish Creek (140500010407) 1978-10-01 2022-06-15
## 6    SNTL    CO    tower  Fish Creek (140500010407) 1978-10-01 2022-06-15
##   latitude longitude elev county site_id       date snow_water_equivalent
## 1    40.54   -106.68 3237  Routt     825 1979-08-13                    NA
## 2    40.54   -106.68 3237  Routt     825 1979-08-14                    NA
## 3    40.54   -106.68 3237  Routt     825 1979-08-15                    NA
## 4    40.54   -106.68 3237  Routt     825 1979-08-16                    NA
## 5    40.54   -106.68 3237  Routt     825 1979-08-17                    NA
## 6    40.54   -106.68 3237  Routt     825 1979-08-18                    NA
##   precipitation_cumulative temperature_max temperature_min temperature_mean
## 1                      841              NA              NA               NA
## 2                      861              NA              NA               NA
## 3                       NA              NA              NA               NA
## 4                       NA              NA              NA               NA
## 5                       NA              NA              NA               NA
## 6                       NA              NA              NA               NA
##   precipitation
## 1            20
## 2            NA
## 3            NA
## 4            NA
## 5            NA
## 6            NA

snotel_825_tower$Date <- as.Date(snotel_825_tower$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

snotel_825_clean <- snotel_825_tower %>% # filter for the timeframe
  filter(Date >= "1978-10-01" & Date <= "2021-09-30") %>%
  filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

adding water day using difftime (SUPER COOL. example from this)

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_825_clean <- snotel_825_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))


write.csv(snotel_825_clean,"C:/Users/13074/Documents/ESS580/thesis_project/thesis_research/data_clean/snotel_825_clean.csv", row.names = FALSE)

Figure check

ggplot(snotel_825_clean, aes(x = Date, y = temperature_mean, filter(waterYear = 1990))) + #this filter didn't work.....
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

#Check for outliers....

#dygraph

temp_xts <- xts(snotel_825_clean$temperature_mean, order.by = snotel_825_clean$Date)

dygraph(temp_xts) %>%
  dyAxis("y", label = "Daily temperature (°C)") 

Detrending Data

#SF figured out the yearly average by water year

#average water year temperature

yearly_wy_aver <- snotel_825_clean %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver <- yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver <- daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver$aver_day_temp))

#str(daily_wy_aver)

# try to show all years as means. 
daily_wy_aver2 <- daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2$date_temp <- signif(daily_wy_aver2$date_temp,3) #reduce the sig figs



ggplot(daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

temp_xts_2 <- xts(daily_wy_aver2$date_temp, order.by = as.Date(daily_wy_aver2$waterDay))

dygraph(temp_xts_2) %>%
  dyAxis("y", label = "(°C)") 

Day of year average temperature for the 1987-2021 period of record for Tower SNOTEL site.

The temperature plateaus in summer

From 266 to 340, there are strange temperature plateaus in the figures.

# try to show all years as means. 
plateau_daily_wy_aver <- daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  filter(waterDay == 266) %>% 
  group_by(waterDay) %>% 
  mutate(date_temp = mean(temperature_mean)) %>% 
  as.tibble()

head(plateau_daily_wy_aver$date_temp, 1)  

## [1] 10.18286

Expanding significant figures fixed it.

All years vs 1987 & 2021

# 87 and 21 only. Not useful
daily_wy_aver3 <- daily_wy_aver2 %>% 
  filter(waterYear == "1987" | waterYear == "2021") %>%
  select(waterDay, waterYear, temperature_mean, daymonth, date_temp) %>%
  group_by(daymonth)

# ggplot(daily_wy_aver3, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
#   geom_line(size= 0.7) +
#   geom_line(aes(y = temperature_mean, color = waterYear)) +
#   theme_few() +
#   #geom_smooth(method = "lm", se=FALSE) +
#   #scale_colour_identity() +
#   scale_color_manual(name = "Water Year", values = c("blue", "red")) 
#   ylab('Daily temperature (°C)') + 
#   xlab('Day of water year')

ggplot(daily_wy_aver3) + #, color = waterYear)) +
  geom_line(aes(x = waterDay, y = date_temp), group = 1, size= 0.7) +
  geom_line(aes(x = waterDay, y = temperature_mean, group = 1, colour = waterYear)) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  #scale_colour_identity() +
  
  #scale_color_manual(name = "Water Year") 
  ylab('Daily temperature (°C)') + 
  xlab('Day of water year')

# ??geom_line  
#   
# scale_color_manual(name = "Colors", 
#                      values = c("a" = "blue", "b" = "red")  
#   
# color_group <- c("blue","black")
#   
# ?scale_discrete_identity()

Day of year average temperature and daily 1987 & 2021 temperatures for the 1987-2021 period of record for Tower SNOTEL site. legend is bad.

Calculating variance

#Calculating residuals, WY 1987 (first full year on record for 825)

residuals_1987 <- daily_wy_aver %>% 
  group_by(waterYear) %>% 
  filter(waterYear == 1987) %>% 
  mutate(residual = (aver_ann_temp-all_ave_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

residuals_2021 <- daily_wy_aver %>% 
  group_by(waterYear) %>% 
  filter(waterYear == 2021) %>% 
  mutate(residual = (aver_ann_temp-all_ave_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

Figures

1987

ggplot(residuals_1987, aes(residual)) +
  geom_histogram(alpha = 0.5, color = "blue", bins = 15, binwidth = 2) +
  theme_base() +
  xlab("Temperature (°C) Above or Below Mean") +
  ylab("Count")

Distribution of daily temperature variance for 1987 compared to the 1987-2021 period of record.

ggplot(residuals_1987, aes(deviation)) +
geom_histogram(alpha = 0.5, color = "blue", bins = 15, binwidth = 1) +
  scale_x_continuous(breaks = seq(0, 15, by=5)) +
  theme_base() +
  xlab("Diurnal Temperature Fluctuation (°C)") +
  ylab("Count")

Distribution of 1987 diurnal temperature fluctuation from the 1987-2021 period of record.

sum(residuals_1987$deviation, na.rm = TRUE)

## [1] 865.0682

865.0682

## [1] 865.0682

residuals_1987 <- residuals_1987 %>% 
  mutate(prcnt = deviation/865.0682*100)

sum(residuals_1987$prcnt, na.rm = TRUE)

## [1] 100

ggplot(residuals_1987, aes(prcnt)) +
geom_histogram(alpha = 0.5, color = "blue", bins = 15, binwidth = 0.1) +
  #scale_x_continuous(breaks = seq(0, 15, by=5)) +
  theme_base() +
  xlab("Diurnal Temperature Fluctuation (°C)") +
  ylab("Percent of Year")

ggplot(residuals_1987) + #, color = waterYear)) +
  geom_jitter(aes(x = deviation, y = prcnt), width = 1, size=1) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  #scale_colour_identity() +
  #scale_color_manual(name = "Water Year") 
  ylab('Percent of Year') + 
  xlab('Deviation (°C)')

Percent of year 1987 for the deviation. This… doesn’t look right.

2021

ggplot(residuals_2021, aes(residual)) +
  geom_histogram(alpha = 0.5, color = "blue", binwidth = 2) +
  theme_base() +
  xlab("Temperature (°C) Above or Below Mean") +
  ylab("Count")

Distribution of daily temperature variance for 2021 compared to the 1987-2021 period of record.

ggplot(residuals_2021, aes(deviation)) +
  geom_histogram(alpha = 0.5, color = "blue", bins = 15, binwidth = 1) +
  scale_x_continuous(breaks = seq(0, 15, by=5)) +
  theme_base() +
  xlab("Diurnal Temperature Fluctuation (°C)") +
  ylab("Count")

# 
# ggplot(residuals_2021, aes(deviation)) +
# geom_histogram(binwidth=1,
#                  center = 6,
#                  aes(col=I("blue"))) +
#   scale_x_continuous(breaks = seq(0, 15, by=5))
#   scale_x_continuous(breaks=seq(1,max(residuals_2021$deviation) + 1, by = 2))

Distribution of 2021 diurnal temperature fluctuation from the 1987-2021 period of record. Scale not the same as 1987

sum(residuals_2021$deviation, na.rm = TRUE)

## [1] 832.3876

# 832.3876

residuals_2021 <- residuals_2021 %>% 
  mutate(prcnt = deviation/832.3876*100)

sum(residuals_2021$prcnt, na.rm = TRUE)

## [1] 100

ggplot(residuals_2021, aes(prcnt)) +
geom_histogram(alpha = 0.5, color = "blue", bins = 15, binwidth = 0.1) +
  #scale_x_continuous(breaks = seq(0, 15, by=5)) +
  theme_base() +
  xlab("Diurnal Temperature Fluctuation (°C)") +
  ylab("Percent of Year")

# This looks weird. Maybe a scatter plot is more appropriate.

ggplot(residuals_2021, aes(prcnt)) +
geom_histogram(alpha = 0.5, color = "blue", bins = 15, binwidth = 0.1) +
  #scale_x_continuous(breaks = seq(0, 15, by=5)) +
  theme_base() +
  xlab("Diurnal Temperature Fluctuation (°C)") +
  ylab("Percent of Year")

ggplot(residuals_2021) + #, color = waterYear)) +
  geom_jitter(aes(x = deviation, y = prcnt), width = 1, size=1) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  #scale_colour_identity() +
  #scale_color_manual(name = "Water Year") 
  ylab('Percent of Year') + 
  xlab('Deviation (°C)')

Percent of year 2021 for the deviaation. This… doesn’t look right.

Mann-Kendall & Sen’s Slope

Mann Kendal for Tower site daily temperatures

# excluding dates where no temperature data were recorded. 

str(daily_wy_aver)

## grouped_df [12,779 x 25] (S3: grouped_df/tbl_df/tbl/data.frame)
##  $ network                 : chr [1:12779] "SNTL" "SNTL" "SNTL" "SNTL" ...
##  $ state                   : chr [1:12779] "CO" "CO" "CO" "CO" ...
##  $ site_name               : chr [1:12779] "tower " "tower " "tower " "tower " ...
##  $ description             : chr [1:12779] "Fish Creek (140500010407)" "Fish Creek (140500010407)" "Fish Creek (140500010407)" "Fish Creek (140500010407)" ...
##  $ start                   : Date[1:12779], format: "1978-10-01" "1978-10-01" ...
##  $ end                     : Date[1:12779], format: "2022-06-15" "2022-06-15" ...
##  $ latitude                : num [1:12779] 40.5 40.5 40.5 40.5 40.5 ...
##  $ longitude               : num [1:12779] -107 -107 -107 -107 -107 ...
##  $ elev                    : num [1:12779] 3237 3237 3237 3237 3237 ...
##  $ county                  : chr [1:12779] "Routt" "Routt" "Routt" "Routt" ...
##  $ site_id                 : num [1:12779] 825 825 825 825 825 825 825 825 825 825 ...
##  $ date                    : chr [1:12779] "1985-09-30" "1985-10-01" "1985-10-02" "1985-10-03" ...
##  $ snow_water_equivalent   : int [1:12779] 15 0 0 0 0 0 0 20 66 91 ...
##  $ precipitation_cumulative: int [1:12779] 1524 0 0 0 0 0 0 20 48 71 ...
##  $ temperature_max         : num [1:12779] 0 0 0 0 0 0 0 0 0 0 ...
##  $ temperature_min         : num [1:12779] 0 0 0 0 0 0 0 0 0 0 ...
##  $ temperature_mean        : num [1:12779] 0 0 0 0 0 0 0 0 0 0 ...
##  $ precipitation           : int [1:12779] 0 0 0 0 0 0 20 28 23 0 ...
##  $ Date                    : Date[1:12779], format: "1985-09-30" "1985-10-01" ...
##  $ waterYear               : num [1:12779] 1985 1986 1986 1986 1986 ...
##  $ daymonth                : chr [1:12779] "30-09" "01-10" "02-10" "03-10" ...
##  $ waterDay                : int [1:12779] 365 1 2 3 4 5 6 7 8 9 ...
##  $ aver_ann_temp           : num [1:12779] 0 1.69 1.69 1.69 1.69 ...
##  $ aver_day_temp           : num [1:12779] 0.362 0.478 0.478 0.479 0.431 ...
##  $ all_ave_temp            : num [1:12779] 0.397 0.397 0.397 0.397 0.397 ...
##  - attr(*, "groups")= tibble [366 x 2] (S3: tbl_df/tbl/data.frame)
##   ..$ daymonth: chr [1:366] "01-01" "01-02" "01-03" "01-04" ...
##   ..$ .rows   : list<int> [1:366] 
##   .. ..$ : int [1:36] 88 443 795 1161 1518 1877 2235 2597 2916 3259 ...
##   .. ..$ : int [1:35] 118 474 826 1192 1549 1908 2266 2628 2946 3290 ...
##   .. ..$ : int [1:33] 144 501 855 1220 1577 1936 2295 2656 2972 3318 ...
##   .. ..$ : int [1:35] 175 532 886 1251 1608 1967 2326 2687 3003 3349 ...
##   .. ..$ : int [1:36] 205 562 916 1281 1638 1997 2356 2717 3032 3379 ...
##   .. ..$ : int [1:36] 236 593 947 1312 1669 2027 2387 2748 3063 3398 ...
##   .. ..$ : int [1:36] 265 623 977 1341 1697 2057 2417 2778 3083 3428 ...
##   .. ..$ : int [1:35] 295 654 1008 1370 1726 2448 2809 3113 3459 3822 ...
##   .. ..$ : int [1:33] 325 684 1039 1399 1756 2113 2840 3144 3489 3853 ...
##   .. ..$ : int [1:34] 2 355 703 1069 1427 1785 2143 2506 2869 3517 ...
##   .. ..$ : int [1:35] 33 383 734 1100 1457 1816 2174 2536 2891 3198 ...
##   .. ..$ : int [1:34] 62 413 764 1130 1487 1846 2204 2566 3228 3578 ...
##   .. ..$ : int [1:36] 89 444 796 1162 1519 1878 2236 2598 2917 3260 ...
##   .. ..$ : int [1:36] 119 475 827 1193 1550 1909 2267 2629 2947 3291 ...
##   .. ..$ : int [1:34] 145 502 856 1221 1578 1937 2296 2657 2973 3319 ...
##   .. ..$ : int [1:35] 176 533 887 1252 1609 1968 2327 2688 3004 3350 ...
##   .. ..$ : int [1:36] 206 563 917 1282 1639 1998 2357 2718 3033 3380 ...
##   .. ..$ : int [1:36] 237 594 948 1313 1670 2028 2388 2749 3064 3399 ...
##   .. ..$ : int [1:36] 266 624 978 1342 1698 2058 2418 2779 3084 3429 ...
##   .. ..$ : int [1:35] 296 655 1009 1371 1727 2449 2810 3114 3460 3823 ...
##   .. ..$ : int [1:35] 326 685 1040 1400 1757 2114 2479 2841 3145 3490 ...
##   .. ..$ : int [1:34] 3 356 704 1070 1428 1786 2144 2507 2870 3518 ...
##   .. ..$ : int [1:35] 34 384 735 1101 1458 1817 2175 2537 2892 3199 ...
##   .. ..$ : int [1:35] 63 414 765 1131 1488 1847 2205 2567 3229 3579 ...
##   .. ..$ : int [1:36] 90 445 797 1163 1520 1879 2237 2599 2918 3261 ...
##   .. ..$ : int [1:36] 120 476 828 1194 1551 1910 2268 2630 2948 3292 ...
##   .. ..$ : int [1:34] 146 503 857 1222 1579 1938 2297 2658 2974 3320 ...
##   .. ..$ : int [1:35] 177 534 888 1253 1610 1969 2328 2689 3005 3351 ...
##   .. ..$ : int [1:35] 207 564 918 1283 1640 2358 2719 3034 3381 3732 ...
##   .. ..$ : int [1:36] 238 595 949 1314 1671 2029 2389 2750 3065 3400 ...
##   .. ..$ : int [1:36] 267 625 979 1343 1699 2059 2419 2780 3085 3430 ...
##   .. ..$ : int [1:36] 297 656 1010 1372 1728 2088 2450 2811 3115 3461 ...
##   .. ..$ : int [1:34] 327 1041 1401 1758 2115 2480 2842 3146 3491 3855 ...
##   .. ..$ : int [1:33] 4 357 705 1071 1429 1787 2145 2508 3519 4205 ...
##   .. ..$ : int [1:35] 35 385 736 1102 1459 1818 2176 2538 2893 3200 ...
##   .. ..$ : int [1:35] 64 415 766 1132 1489 1848 2206 2568 3230 3580 ...
##   .. ..$ : int [1:36] 91 446 798 1164 1521 1880 2238 2600 2919 3262 ...
##   .. ..$ : int [1:36] 121 477 829 1195 1552 1911 2269 2631 2949 3293 ...
##   .. ..$ : int [1:33] 147 504 858 1223 1580 1939 2298 2659 2975 3321 ...
##   .. ..$ : int [1:36] 178 535 889 1254 1611 1970 2329 2690 3006 3352 ...
##   .. ..$ : int [1:36] 208 565 919 1284 1641 1999 2359 2720 3035 3382 ...
##   .. ..$ : int [1:36] 239 596 950 1315 1672 2030 2390 2751 3066 3401 ...
##   .. ..$ : int [1:35] 268 626 980 1344 2060 2420 2781 3086 3431 3794 ...
##   .. ..$ : int [1:36] 298 657 1011 1373 1729 2089 2451 2812 3116 3462 ...
##   .. ..$ : int [1:35] 328 686 1042 1402 1759 2116 2481 2843 3147 3492 ...
##   .. ..$ : int [1:35] 5 358 706 1072 1430 1788 2146 2509 2871 3170 ...
##   .. ..$ : int [1:36] 36 386 737 1103 1460 1819 2177 2539 2894 3201 ...
##   .. ..$ : int [1:35] 65 416 767 1133 1490 1849 2207 2569 3231 3581 ...
##   .. ..$ : int [1:34] 92 447 799 1165 1522 1881 2239 2601 2920 3263 ...
##   .. ..$ : int [1:35] 122 478 830 1196 1553 1912 2270 2632 2950 3294 ...
##   .. ..$ : int [1:34] 148 505 859 1224 1581 1940 2299 2660 2976 3322 ...
##   .. ..$ : int [1:36] 179 536 890 1255 1612 1971 2330 2691 3007 3353 ...
##   .. ..$ : int [1:36] 209 566 920 1285 1642 2000 2360 2721 3036 3383 ...
##   .. ..$ : int [1:36] 240 597 951 1316 1673 2031 2391 2752 3067 3402 ...
##   .. ..$ : int [1:35] 269 627 981 1345 2061 2421 2782 3087 3432 3795 ...
##   .. ..$ : int [1:35] 299 658 1012 1730 2090 2452 2813 3117 3463 3826 ...
##   .. ..$ : int [1:35] 329 687 1043 1403 1760 2117 2482 2844 3148 3493 ...
##   .. ..$ : int [1:35] 6 359 707 1073 1431 1789 2147 2510 2872 3171 ...
##   .. ..$ : int [1:35] 37 387 738 1104 1461 1820 2178 2540 3202 3552 ...
##   .. ..$ : int [1:35] 66 417 768 1134 1491 1850 2208 2570 3232 3582 ...
##   .. ..$ : int [1:33] 93 448 800 1166 1523 1882 2240 2602 2921 3264 ...
##   .. ..$ : int [1:36] 123 479 831 1197 1554 1913 2271 2633 2951 3295 ...
##   .. ..$ : int [1:33] 149 506 860 1225 1582 1941 2300 2661 2977 3323 ...
##   .. ..$ : int [1:35] 180 537 891 1256 1613 1972 2331 2692 3008 3354 ...
##   .. ..$ : int [1:36] 210 567 921 1286 1643 2001 2361 2722 3037 3384 ...
##   .. ..$ : int [1:35] 241 598 952 1674 2032 2392 2753 3068 3403 3766 ...
##   .. ..$ : int [1:35] 270 628 982 1346 1700 2062 2422 2783 3088 3433 ...
##   .. ..$ : int [1:36] 300 659 1013 1374 1731 2091 2453 2814 3118 3464 ...
##   .. ..$ : int [1:34] 330 688 1044 1761 2118 2483 2845 3149 3494 3858 ...
##   .. ..$ : int [1:33] 7 360 708 1074 1432 1790 2148 2511 3172 3522 ...
##   .. ..$ : int [1:35] 38 388 739 1105 1462 1821 2179 2541 3203 3553 ...
##   .. ..$ : int [1:34] 418 769 1135 1492 1851 2209 2571 3233 3583 3909 ...
##   .. ..$ : int [1:35] 94 449 801 1167 1524 1883 2241 2603 2922 3265 ...
##   .. ..$ : int [1:36] 124 480 832 1198 1555 1914 2272 2634 2952 3296 ...
##   .. ..$ : int [1:33] 150 507 861 1226 1583 1942 2301 2662 2978 3324 ...
##   .. ..$ : int [1:35] 181 538 892 1257 1614 1973 2332 2693 3009 3355 ...
##   .. ..$ : int [1:36] 211 568 922 1287 1644 2002 2362 2723 3038 3385 ...
##   .. ..$ : int [1:35] 242 599 953 1317 1675 2033 2393 2754 3404 3767 ...
##   .. ..$ : int [1:36] 271 629 983 1347 1701 2063 2423 2784 3089 3434 ...
##   .. ..$ : int [1:36] 301 660 1014 1375 1732 2092 2454 2815 3119 3465 ...
##   .. ..$ : int [1:34] 331 1045 1404 1762 2119 2484 2846 3150 3495 3859 ...
##   .. ..$ : int [1:33] 8 361 709 1075 1433 1791 2149 2512 3173 3523 ...
##   .. ..$ : int [1:35] 39 389 740 1106 1463 1822 2180 2542 3204 3554 ...
##   .. ..$ : int [1:35] 419 770 1136 1493 1852 2210 2572 2895 3234 3584 ...
##   .. ..$ : int [1:35] 95 450 802 1168 1525 1884 2242 2604 2923 3266 ...
##   .. ..$ : int [1:35] 125 481 833 1199 1556 1915 2273 2635 2953 3297 ...
##   .. ..$ : int [1:32] 151 508 862 1227 1584 1943 2302 2663 2979 3325 ...
##   .. ..$ : int [1:36] 182 539 893 1258 1615 1974 2333 2694 3010 3356 ...
##   .. ..$ : int [1:36] 212 569 923 1288 1645 2003 2363 2724 3039 3386 ...
##   .. ..$ : int [1:35] 243 600 954 1318 1676 2034 2394 2755 3405 3768 ...
##   .. ..$ : int [1:34] 272 630 984 1348 1702 2064 2424 2785 3435 3798 ...
##   .. ..$ : int [1:36] 302 661 1015 1376 1733 2093 2455 2816 3120 3466 ...
##   .. ..$ : int [1:34] 332 689 1046 1405 1763 2120 2485 2847 3151 3496 ...
##   .. ..$ : int [1:33] 9 362 710 1076 1434 1792 2150 2513 3174 3524 ...
##   .. ..$ : int [1:35] 40 390 741 1107 1464 1823 2181 2543 3205 3555 ...
##   .. ..$ : int [1:35] 420 771 1137 1494 1853 2211 2573 2896 3235 3585 ...
##   .. ..$ : int [1:34] 96 451 803 1169 1526 1885 2243 2605 3267 3617 ...
##   .. ..$ : int [1:35] 482 834 1200 1557 1916 2274 2636 2954 3298 3648 ...
##   .. ..$ : int [1:33] 152 509 863 1228 1585 1944 2303 2664 2980 3326 ...
##   .. .. [list output truncated]
##   .. ..@ ptype: int(0) 
##   ..- attr(*, ".drop")= logi TRUE

daily_wy_aver4 <- daily_wy_aver2 %>% 
  filter(Date > "1986-07-17")

tower_mk <- mk.test(daily_wy_aver4$temperature_mean)
print(tower_mk)

## 
##  Mann-Kendall trend test
## 
## data:  daily_wy_aver4$temperature_mean
## z = 10.156, n = 12499, p-value < 2.2e-16
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 4.731005e+06 2.169853e+11 6.066612e-02

Sen’s Slope for Tower site daily temperatures

tower_sens <- sens.slope(daily_wy_aver4$temperature_mean)

print(tower_sens)

## 
##  Sen's slope
## 
## data:  daily_wy_aver4$temperature_mean
## z = 10.156, n = 12499, p-value < 2.2e-16
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.0001879699 0.0002786291
## sample estimates:
##  Sen's slope 
## 0.0002328967

SNOTEL Site corrections

“The Morrissey method also shows more variations in temperature trends between sites than the Oyler method, which dampens site-to-site variability in temperature trends.

Next steps / Further questions

SNOTEL

Need to investigate site specific issues related to disturbances (general SNOTEL issues). Unsure of methods

Having issues with day-of-year and water year interactions. This keeps me from showing mean trendlines over the whole period of record as R does not want to recognize the x axis as anything other than numerical. Resolved

Would using a pivot table be better when displaying temperature for 240 sites? (What data management is best for that much data?) Doesn’t look like it’s needed currently

How to include oscillation influence into the data?

Beyond SNOTEL

What other data / site types do I want to include?

Coding

I am sure Matt would tell me to automate this. Not quite sure how, but it seems likely I could find a list of SNOTEL sites and using spatial data, exclude sites not within the CRB boundary. From there, I could write a for loop & function to use the snotelr function to download data. Jack Reuland helped with this idea below. Now what to do with the large dataframe….

All stations:

Per Jack’s recommendations, BUT, while there are 240 stations, I could only find 166 while clipping the NRCS SNOTEL points and a CRB polygon from CAPGISadmin. The NRCS stated the sites given were active, so perhaps there are 74 inactive sites?

TableOfSNOTEL_StudyArea <- read.csv("C:/Users/13074/Desktop/MS_Watershed_Science/Literature/SNOTEL_sites_1_2.csv", header = TRUE)

studyarea_site_id <- as.vector(TableOfSNOTEL_StudyArea$site_id)

SNOTEL_StudyArea <- snotel_download(site_id = c(studyarea_site_id), internal = TRUE)

## Error in wdman::phantomjs(verbose = FALSE) : 
##   PhantomJS signals port = 4567 is already in use.

write.csv(SNOTEL_StudyArea,"C:/Users/13074/Documents/ESS580/thesis_project/thesis_research/data_raw/SNOTEL_StudyArea.csv", row.names = FALSE) #write in the raw data

head(SNOTEL_StudyArea)

##   network state       site_name                             description
## 1    SNTL    CO elkhead divide  Headwaters Elkhead Creek (140500010601)
## 2    SNTL    CO elkhead divide  Headwaters Elkhead Creek (140500010601)
## 3    SNTL    CO elkhead divide  Headwaters Elkhead Creek (140500010601)
## 4    SNTL    CO elkhead divide  Headwaters Elkhead Creek (140500010601)
## 5    SNTL    CO elkhead divide  Headwaters Elkhead Creek (140500010601)
## 6    SNTL    CO elkhead divide  Headwaters Elkhead Creek (140500010601)
##        start        end latitude longitude elev county site_id       date
## 1 2016-09-01 2022-06-21     40.8    -107.1 2682  Routt    1252 2016-09-01
## 2 2016-09-01 2022-06-21     40.8    -107.1 2682  Routt    1252 2016-09-02
## 3 2016-09-01 2022-06-21     40.8    -107.1 2682  Routt    1252 2016-09-03
## 4 2016-09-01 2022-06-21     40.8    -107.1 2682  Routt    1252 2016-09-04
## 5 2016-09-01 2022-06-21     40.8    -107.1 2682  Routt    1252 2016-09-05
## 6 2016-09-01 2022-06-21     40.8    -107.1 2682  Routt    1252 2016-09-06
##   snow_water_equivalent precipitation_cumulative temperature_max
## 1                    NA                       NA            23.6
## 2                     0                      795            20.6
## 3                     0                      798            17.8
## 4                     0                      823            15.8
## 5                     0                      833            19.0
## 6                     0                      836            19.6
##   temperature_min temperature_mean precipitation
## 1            11.0             16.4            NA
## 2             9.8             14.1             3
## 3             8.1             11.3            25
## 4             6.8             10.9            10
## 5             6.4             12.3             3
## 6             6.4             12.3             0

Coop Data

NOAA/NWS Cooperative Observer Network

Two stations near steamboat, according to this site. 057936 02 SBTC2 UNITED STATES CO ROUTT +7 STEAMBOAT SPRINGS 40 30 17 -106 51 58 6636 057942 02 SSPC2 UNITED STATES CO ROUTT +7 STEAMBOAT SPRINGS 1 W 40 29 00 -106 51 00 6700

Downloading from the Iowa Environmental Mesonet. NOAA package not working, maybe due to user error.

# coops_search(station_name = 057936,product = "air_temperature") # Handy R package did not work.

coop_steamboat <- read.csv("C:/Users/13074/Documents/ESS580/thesis_project/thesis_research/data_raw/steamboat_nwscoop_1900_2022.csv", header = TRUE)

coop_steamboat$Date <- mdy(coop_steamboat$day) #change date from character to date format, capitalize to work with Water year functon from NWIS.

coop_sb_cleanish <- coop_steamboat %>% 
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#str(coop_steamboat)

sb_coop_temp_xts <- xts(coop_sb_cleanish$avg_T_c, order.by = coop_sb_cleanish$Date)

dygraph(sb_coop_temp_xts) %>%
  dyAxis("y", label = "Daily temperature (°C)") 

Notes from SF, 22/6:

Some initial thoughts: 1) I assume that frequency is the raw count, i.e., number of days? At some point, a % of the year may be better; 2) COOP data, such as at Steamboat, should have a much longer period of record if you want to examine change over time - setting up the analysis is a good start; 3) there is an inconsistency with the SNOTEL data (see https://doi.org/10.1029/2019WR025921); 4) we (https://doi.org/10.1029/2002WR001512, etc.) used 240 stations (metadata at https://doi.org/10.1029/2009WR007835), and I know that are a few are inactive, but likely not 74 - some are outside of the basin (https://agupubs.onlinelibrary.wiley.com/cms/asset/2af0fc20-a948-455d-bf0e-8f6e7dcc7ad9/wrcr12217-fig-0001.png) but still representative (a few can be removed, in hindsight, like the ones in the Laramie Range).

Notes, Monday 27/6

Figured out the water day issue and applied it to the Day of water year figures (this was an issue that I spent waaaaaay too much time on). Went the wrong way with pivot tables, day of year, and using average annual temp for the mean aily figure. Going to use Excel to check if correct. Goal for Tuesday: Determine daily percentages.

Notes, Thursday 30/6

Need to discuss methods of identifying temperature shifts due to sensor changes.

Not sure how to show percent of days for data.

Corrected Variance 19/7

The above chunks of code may be incorrect; instead of subtracting the mean temperature for the period of record from the mean annual temperature for each water year, I should switch that and subtract the mean annual temperature for each water year from mean temperature for the period of record to find the residuals:

(all_ave_temp - aver_ann_temp) + temperature_mean - aver_day_temp The average of the residuals should be zero.

residuals_1987_2 <- daily_wy_aver %>% 
  group_by(waterYear) %>% 
  filter(waterYear == 1987) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(residuals_1987_2$residual)

## [1] 0.001015426

residuals_2021_2 <- daily_wy_aver %>% 
  group_by(waterYear) %>% 
  filter(waterYear == 2021) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(residuals_2021_2$residual)

## [1] -0.0007127652

Based on the calculated means, this looks correct.

1987

ggplot(residuals_1987_2, aes(residual)) +
  geom_histogram(alpha = 0.5, color = "blue", bins = 15, binwidth = 2) +
  theme_base() +
  xlab("Temperature (°C) Above or Below Mean") +
  ylab("Count")

Distribution of daily temperature variance for 1987 compared to the 1987-2021 period of record.

ggplot(residuals_1987_2, aes(deviation)) +
geom_histogram(alpha = 0.5, color = "blue", bins = 15, binwidth = 1) +
  scale_x_continuous(breaks = seq(0, 15, by=5)) +
  theme_base() +
  xlab("Diurnal Temperature Fluctuation (°C)") +
  ylab("Count")

Distribution of 1987 diurnal temperature fluctuation from the 1987-2021 period of record.

2021

ggplot(residuals_2021_2, aes(residual)) +
  geom_histogram(alpha = 0.5, color = "blue", binwidth = 2) +
  theme_base() +
  xlab("Temperature (°C) Above or Below Mean") +
  ylab("Count")

Distribution of daily temperature variance for 2021 compared to the 1987-2021 period of record.

ggplot(residuals_2021_2, aes(deviation)) +
  geom_histogram(alpha = 0.5, color = "blue", bins = 15, binwidth = 1) +
  scale_x_continuous(breaks = seq(0, 15, by=5)) +
  theme_base() +
  xlab("Diurnal Temperature Fluctuation (°C)") +
  ylab("Count")

# 
# ggplot(residuals_2021, aes(deviation)) +
# geom_histogram(binwidth=1,
#                  center = 6,
#                  aes(col=I("blue"))) +
#   scale_x_continuous(breaks = seq(0, 15, by=5))
#   scale_x_continuous(breaks=seq(1,max(residuals_2021$deviation) + 1, by = 2))

Distribution of 2021 diurnal temperature fluctuation from the 1987-2021 period of record. Scale not the same as 1987

Standard Deviation

To figure out the standard deviation for each year, I want the “residual” for each daily value.

The standard deviation will be the daily residual minus the mean of the residuals by water year, summed and squared, then divided by the number of observations minus one. The square root of the resulting value of which is thus the standard deviation for the water year.

standard_dev <- daily_wy_aver %>% 
  group_by(waterYear) %>% 
  filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_87 <- standard_dev %>% 
  filter(waterYear == 1987) %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
           mutate(sd_1 = residual-resid_mean)

standard_dev_87 <- standard_dev_87 %>%
  group_by(waterYear) %>%
  mutate(sd_2 = (((sum((sd_1)^2))/((sum(tabulate(standard_dev_87$waterDay))))))^(0.5))

standard_dev_88 <- standard_dev %>% 
  filter(waterYear == 1988) %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
           mutate(sd_1 = residual-resid_mean)

standard_dev_88 <- standard_dev_88 %>%
  group_by(waterYear) %>%
  mutate(sd_2 = (((sum((sd_1)^2))/((sum(tabulate(standard_dev_87$waterDay))))))^(0.5))

standard_dev_89 <- standard_dev %>% 
  filter(waterYear == 198) %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
           mutate(sd_1 = residual-resid_mean)

standard_dev_89 <- standard_dev_89 %>%
  group_by(waterYear) %>%
  mutate(sd_2 = (((sum((sd_1)^2))/((sum(tabulate(standard_dev_87$waterDay))))))^(0.5))

Standard Deviation each year (likely not correct)

standard_dev_all <- standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all <- standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')

waterYear	sd_2
1987	8.869131
1988	10.141239
1989	9.503716
1990	9.373486
1991	9.479342
1992	8.453773
1993	8.671267
1994	9.749269
1995	8.844523
1996	9.223301
1997	9.523966
1998	9.430494
1999	8.883579
2000	8.937816
2001	9.612584
2002	9.783074
2003	9.303952
2004	8.779754
2005	8.155340
2006	9.030542
2007	9.099075
2008	9.487077
2009	8.307930
2010	8.948885
2011	8.857272
2012	8.984175
2013	9.367552
2014	8.635253
2015	8.092495
2016	8.777249
2017	8.579684
2018	8.620667
2019	8.849289
2020	9.314919
2021	9.282680

ggplot(standard_dev_all, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')